home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FishMarket 1.0
/
FishMarket v1.0.iso
/
fishies
/
301-325
/
disk_321
/
planets
/
moon_pos.c
< prev
next >
Wrap
C/C++ Source or Header
|
1992-05-06
|
8KB
|
237 lines
/*
* moon-pos.c - Determine the position of the moon.
*
* Author: Jim Cobb
* jcobb@utah.edu.cs
* Computer Science Dept.
* University of Utah
* Date: Fri Mar 25 1988
*
* Copyright (c) 1988, Jim Cobb
* All rights reserved
*
* (This copyright prevents you from selling the program - the
* author grants anyone the right to copy and install the program
* on any machine it will run on)
*/
/* Reference: Jean Meeus, " Astronomical Formulae for Calculators,"
* chapter 30. Meeus states that this procedure has an accuracy of 10"
* in the longitude of the moon, 3" in its latitude and 0.2" in its
* parallax.
*/
#include <stdio.h>
#include <math.h>
/* Trigonometric functions in degrees. */
#ifndef PI
# define PI 3.1415926535897932385 /* Should be in math.h . */
#endif /* PI */
#define RAD (PI/180.0)
#define sind( x ) ( sin( RAD * (x) ) )
#define cosd( x ) ( cos( RAD * (x) ) )
#define MAGMOON -9.9
extern double range();
/* Moon_pos takes a time argument, t, the number of Julian centuries
* from 1900 January 0.5 ET, and the obliquity of the ecliptic, epli.
* It adds information about the moon's position to the logfile.
*/
moon_pos( t, epli )
double t, epli;
{
double t_2 = t*t,
t_3 = t*t*t,
e, e_2,
l_prime, /* Moon's longitude. */
m, /* Sun's anomaly. */
m_prime, /* Moon's anomaly. */
d, /* Moon's elongation. */
f, /* Distance of Moon from its ascending node. */
omega, /* Longitude of Moon's ascending node. */
omega1, omega2,
venus_term, /* Great Venus term. */
lambda, /* Ecliptical longitude of the Moon's center. */
b, beta, /* Ecliptical latitude of the Moon's center. */
parallax, /* Equatorial horizontal parallax of the Moon. */
tmp;
double phase; /* phase of the moon */
static char namestring[] = "Moon, 180.";
/* Compute the mean values for the terms. */
l_prime = 270.434164 + 481267.8831 * t - 0.001133 * t_2 + 0.0000019 * t_3;
m = 358.475833 + 35999.0498 * t - 0.000150 * t_2 - 0.0000033 * t_3;
m_prime = 296.104608 + 477198.8491 * t + 0.009192 * t_2 + 0.0000144 * t_3;
d = 350.737486 + 445267.1142 * t - 0.001436 * t_2 + 0.0000019 * t_3;
f = 11.250889 + 483202.0251 * t - 0.003211 * t_2 - 0.0000003 * t_3;
omega = 259.183275 - 1934.1420 * t + 0.002078 * t_2 + 0.0000022 * t_3;
/* Additive terms. */
tmp = sind( 51.2 + 20.2 * t );
l_prime += 0.000233 * tmp;
m += -0.001778 * tmp;
m_prime += 0.000817 * tmp;
d += 0.002011 * tmp;
venus_term = 0.003964 *
sind( 346.560 + 132.870 * t - 0.0091731 * t_2 );
l_prime += venus_term;
m_prime += venus_term;
d += venus_term;
f += venus_term;
tmp = sind( omega );
l_prime += 0.001964 * tmp;
m_prime += 0.002541 * tmp;
d += 0.001964 * tmp;
f += -0.024691 * tmp;
f += -0.004328 * sind( omega + 275.05 - 2.30 * t );
e = 1 - 0.002495 * t - 0.00000752 * t_2;
e_2 = e*e;
/* Bring these angles within 0 to 360 degrees. */
m = range( m );
m_prime = range( m_prime );
d = range( d );
f = range( f );
omega = range( omega );
/* broke into pieces so a micro compute compiler could handel it -- Bob L */
/* Calculate lambda, the ecliptical longitude of the Moon's center. */
lambda = l_prime + 6.288750 * sind( m_prime )
+ 1.274018 * sind( 2*d - m_prime )
+ 0.658309 * sind( 2*d )
+ 0.213616 * sind( 2 * m_prime )
- e * 0.185596 * sind( m )
- 0.114336 * sind( 2*f )
+ 0.058793 * sind( 2*d - 2*m_prime )
+ e * 0.057212 * sind( 2*d - m - m_prime );
lambda += 0.053320 * sind( 2*d + m_prime )
+ e * 0.045874 * sind( 2*d - m )
+ e * 0.041024 * sind( m_prime - m )
- 0.034718 * sind( d )
- e * 0.030465 * sind( m + m_prime )
+ 0.015326 * sind( 2*d - 2*f )
- 0.012528 * sind( 2*f + m_prime )
- 0.010980 * sind( 2*f - m_prime );
lambda += 0.010674 * sind( 4*d - m_prime )
+ 0.010034 * sind( 3*m_prime )
+ 0.008548 * sind( 4*d - 2*m_prime )
- e * 0.007910 * sind( m - m_prime + 2*d )
- e * 0.006783 * sind( 2*d + m )
+ 0.005162 * sind( m_prime - d )
+ e * 0.005000 * sind( m + d )
+ e * 0.004049 * sind( m_prime - m + 2*d )
+ 0.003996 * sind( 2*m_prime + 2*d )
+ 0.003862 * sind( 4*d );
lambda += 0.003665 * sind( 2*d - 3*m_prime )
+ e * 0.002695 * sind( 2*m_prime - m )
+ 0.002602 * sind( m_prime - 2*f - 2*d )
+ e * 0.002396 * sind( 2*d - m - 2*m_prime )
- 0.002349 * sind( m_prime + d )
+ e_2 * 0.002249 * sind( 2*d - 2*m )
- e * 0.002125 * sind( 2*m_prime + m )
- e_2 * 0.002079 * sind( 2*m )
+ e_2 * 0.002059 * sind( 2*d - m_prime - 2*m )
- 0.001773 * sind( m_prime + 2*d - 2*f )
- 0.001595 * sind( 2*f + 2*d );
lambda += e * 0.001220 * sind( 4*d - m - m_prime )
- 0.001110 * sind( 2*m_prime + 2*f )
+ 0.000892 * sind( m_prime - 3*d )
- e * 0.000811 * sind( m + m_prime + 2*d )
+ e * 0.000761 * sind( 4*d - m - 2*m_prime )
+ e_2 * 0.000717 * sind( m_prime - 2*m )
+ e_2 * 0.000704 * sind( m_prime - 2*m -2*d )
+ e * 0.000693 * sind( m - 2*m_prime + 2*d )
+ e * 0.000598 * sind( 2*d - m - 2*f )
+ 0.000550 * sind( m_prime + 4*d )
+ 0.000538 * sind( 4*m_prime )
+ e * 0.000521 * sind( 4*d - m )
+ 0.000486 * sind( 2*m_prime - d );
lambda = range( lambda );
/* Calculate beta, the ecliptical latitude of the Moon's center. */
b = 5.128189 * sind( f )
+ 0.280606 * sind( m_prime + f )
+ 0.277693 * sind( m_prime - f )
+ 0.173238 * sind( 2*d - f )
+ 0.055413 * sind( 2*d + f - m_prime )
+ 0.046272 * sind( 2*d - f - m_prime )
+ 0.032573 * sind( 2*d + f )
+ 0.017198 * sind( 2*m_prime + f )
+ 0.009267 * sind( 2*d + m_prime - f )
+ 0.008823 * sind( 2*m_prime - f )
+ e * 0.008247 * sind( 2*d - m - f )
+ 0.004323 * sind( 2*d - f - 2*m_prime );
b += 0.004200 * sind( 2*d + f + m_prime )
+ e * 0.003372 * sind( f - m - 2*d )
+ e * 0.002472 * sind( 2*d + f - m - m_prime )
+ e * 0.002222 * sind( 2*d + f - m )
+ e * 0.002072 * sind( 2*d - f - m - m_prime )
+ e * 0.001877 * sind( f - m + m_prime )
+ 0.001828 * sind( 4*d - f - m_prime )
- e * 0.001803 * sind( f + m )
- 0.001750 * sind( 3*f )
+ e * 0.001570 * sind( m_prime - m - f )
- 0.001487 * sind( f + d )
- e * 0.001481 * sind( f + m + m_prime );
b += e* 0.001417 * sind( f - m - m_prime )
+ e * 0.001350 * sind( f - m )
+ 0.001330 * sind( f - d )
+ 0.001106 * sind( f + 3*m_prime )
+ 0.001020 * sind( 4*d - f )
+ 0.000833 * sind( f + 4*d - m_prime )
+ 0.000781 * sind( m_prime - 3*f )
+ 0.000670 * sind( f + 4*d - 2*m_prime )
+ 0.000606 * sind( 2*d - 3*f )
+ 0.000597 * sind( 2*d + 2*m_prime - f )
+ e * 0.000492 * sind( 2*d + m_prime - m - f );
b += 0.000450 * sind( 2*m_prime - f - 2*d )
+ 0.000439 * sind( 3*m_prime - f )
+ 0.000423 * sind( f + 2*d + 2*m_prime )
+ 0.000422 * sind( 2*d - f - 3*m_prime )
- e * 0.000367 * sind( m + f + 2*d - m_prime )
- e * 0.000353 * sind( m + f + 2*d )
+ 0.000331 * sind( f + 4*d )
+ e * 0.000317 * sind( 2*d + f - m + m_prime )
+ e_2 * 0.000306 * sind( 2*d - 2*m - f )
- 0.000283 * sind( m_prime + 3*f );
omega1 = 0.0004664 * cosd( omega );
omega2 = 0.0000754 * cosd( omega + 275.05 - 2.30 * t );
beta = b * ( 1 - omega1 - omega2 );
/* roughly calculate the phase of the moon.
Added by ccount
Sun longi is about 36000 * T0 + 279
Phase is 0 for full moon */
tmp = 36000.0 * t +279 - lambda;
phase = 180.0 - (tmp - ((double)((int)(tmp/360.0))*360.0));
phase = (phase < 0.0) ? 360 + phase:phase;
sprintf(namestring, "Moon, %3.0f", phase);
/* Transform to right ascension and declination, and output the
* result.
*/
geo_trans( lambda, beta, epli, 0.0, MAGMOON, "PL", namestring );
}